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Abstract 

Our probabilistic analysis sheds light to the following questions: Why do random poly- 
nomials seem to have few, and well separated real roots, on the average? Why do exact 
f^ ' algorithms for real root isolation may perform comparatively well or even better than nu- 

CN ', merical ones? 

^.. ■ We exploit results by Kac, and by Edelman and Kostlan in order to estimate the real 

C^ , root separation of degree d polynomials with i.i.d. coefficients that follow two zero-mean 

normal distributions: for SO (2) polynomials, the i-th coefficient has variance (^), whereas 
for Weyl polynomials its variance is 1/i!. By applying results from statistical physics, we 
obtain the expected (bit) complexity of STURM solver, OBlfd^T), where r is the number of 
real roots and t the maximum coefficient bitsize. Our bounds are two orders of magnitude 
tighter than the record worst case ones. We also derive an output-sensitive bound in the 
worst case. 
^'J ' The second part of the paper shows that the expected number of real roots of a degree 

C/3 , d polynomial in the Bernstein basis is \/2d± 0(1), when the coefficients are i.i.d. variables 

with moderate standard deviation. Our paper concludes with experimental results which 
corroborate our analysis. 

(N 

*^ , Categories and Subject Descriptors: F.2 [Theory of Computation] : Analysis of Algorithms 

f~^ I and Problem Complexity; I.l [Computing Methodology]: Symbolic and algebraic manipulation: 

C^ ' Algorithms 

Keywords: Random polynomial, real-root isolation, Bernstein polynomial, expected com- 
plexity, separation bound 
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1 Introduction 



^ . One of the most important procedures in computer algebra and algebraic algorithms is root 

H I isolation of univariate polynomials. The goal is to compute intervals in the real case, or squares 

in the complex case, that isolate the roots of the polynomial and to compute one such interval, 
or square, for every root. 

We restrict ourselves to exact algorithms, i.e. algorithms that perform arithmetic with ratio- 
nal numbers of arbitrary size. The best known algorithms are subdivision algorithms, based on 
Sturm sequences (sturm), or on Descartes' rule of sign (descartes), or on Descartes' rule and 
the Bernstein basis representation (bernstein). Subdivision algorithms mimic binary search 
and their complexity depends on separation bounds. They are given an initial interval, or com- 
pute one containing all real roots. Then, they repeatedly subdivide it until it is certified that 
zero or one real root is contained in the tested interval. 
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Thanks to important recent progress [3 El [TUl [TI] , the complexity of STURM, DESCARTES and 
BERNSTEIN is, in the worst case, O^id "C ]^ where d is the degree of the polynomial and T the 
maximum coefficient bitsize. The bound holds even when the polynomial is non-squarefree, and 
we also compute (all) the multiplicities. This requires a preprocessing of complexity Oeld t), 
in order to compute the square-free factorization. The new polynomial has coefhcients of size 
0{d + t). The complexity of this stage, although significant in practice, is asymptotically 
dominated. In this paper we consider the behavior of STURM on random polynomials of various 
forms. Our results can be extended to DESCARTES and BERNSTEIN. 

Another important exact solver (cf) is based on the continued fractions expansion of the 
real roots e.g. [Tl[33l[35]. Several variants of this solver exist, depending on the method used to 
compute the partial quotients of the real roots. Assuming the Gauss-Kuzmin distribution holds 
for the real algebraic numbers, it was proven |35) . that the expected complexity is Ogfd'^T ). 
By spreading the roots, the expected complexity becomes O^id "v) |35| . The currently known 
worst-case bound is O^id^l ) [25| . This paper reduces the gap between STURM CF. 

Numerical algorithms compute an approximation, up to a desired accuracy, of all complex 
roots. They can be turned into isolation algorithms by requiring the accuracy to be equal to 
the theoretical worst-case separation bound. The current record is O^id r) and is achieved by 
recursively splitting the polynomial until one obtains linear factors that approximate sufficiently 
the roots |321I27| . It seems that the bounds could be improved to O^id r] with a more sophisti- 
cated splitting process. We should mention that optimal numerical algorithms are very difficult 
to implement. 

Even though the complexity bounds of the exact algorithms are worse than those of the 
numerical ones, recent implementations of the former tend to be competitive, if not superior, in 
practice, e.g. |19l [30l[Tni35| . Our work attempts to provide an explanation for this. There is a 
huge amount of work concerning root isolation and the references stated represent only the tip 
of the iceberg; we encourage the reader to refer to the references. 

Most of the work on random polynomials, which typically concerns polynomials in the mono- 
mial basis, focuses on the number of real roots. Kac's [20j celebrated result estimated the ex- 
pected number of real roots of random polynomials (named after himself) as - log d -|- 0(1), 
when the coefficients are standard normals i.i.d. or uniformly distributed, and d is the degree of 
the polynomial. We refer the reader to e.g. j5l (U |T2] for a historical perspective and to f3] for 
various references. A geometric interpretation of this result and many generalizations appear 
in [9]. We mainly examine S0(2) polynomials, where the l-th coefficient is an i.i.d. Gaussian 
random variable of zero mean and variance (^). According to [S], they are "the most natural 
definition of random polynomials", see also [3¥j . Their expected number of real roots is yd. For 
Weyl polynomials, the l-th coefficient is an i.i.d. Gaussian random variable of zero mean and 
variance 1/1!, and the expected number of real roots is about -vd+ 0(1) where higher-order 
terms are not known to date [;31j. For results on complex roots we refer to e.g. [141113). 

Our first contribution concerns the expected bit complexity of STURM, when the input is 
random polynomials with i.i.d. coefficients; notice that their roots are not independently dis- 
tributed! In other words, we have to go beyond the theory of Kac, and Edelman and Kostlan, 
in order to study the statistical behavior of root differences and, more precisely, the minimum 
absolute difference. We examine SO (2) and Weyl random polynomials, and exploit the relevant 
progress achieved in statistical physics. In fact, these polynomial classes are of particular interest 
in statistical physics because they model zero-crossings in diffusion equations and, eventually, a 
chaotic spin wave- function [H [T^ 131) . The key observation is that, by applying these results, 
we can quantify the correlation between the roots, which is sufficiently weak, but does exist. 
For both classes of polynomials we prove an expected case bit complexity bound of OB(Td t), 
where r is the number of real roots. A close related bound was speculated in [18| . based on 



experimental evidence. 

Our bounds are tighter than those of the worst case by two factors. In the course of this 
analysis, STURM is shown to be output-sensitive, with complexity proportional to the number 
of real roots in the given interval, even in the worst case. A similar bound appeared in |15| . 

Besides polynomials in the monomial basis, polynomials in the Bernstein basis are im- 
portant in many applications, e.g. CAGD and geometric modeling. They are of the form 
^1=0 Q-i(t)''Hl —x)'^^^. For the random polynomials that we consider, at are standard normals 
i.i.d. random variables, that is Gaussians with zero mean and variance one. Such polynomials 
are also important in Brownian motion [21]. In [2], they examine random polynomial systems; 
they also estimate the expected number of real roots of a polynomial in the Bernstein basis as 
yd, when the variance is (^). This left open the case, see also |21| . of smaller variance, that is 
polynomial and not exponential in d. 

Our second contribution is to examine random polynomials in the Bernstein basis of de- 
gree d, with i.i.d. coefficients with mean zero and "moderate" variance 0(1/Y^d/(i(d — i))), for 
d > I > 0. Indeed, we have 1 > y^d/(i(d — I)) > 2/v7td. We prove that the expected number 
of real roots of these polynomials is v2d it 0(1 ). We conclude with experimental results which 
corroborate our analysis, and shows that these polynomials behave like polynomial with vari- 
ance 1. This is the first step towards bounding the expected complexity of solving polynomials 
in the Bernstein basis. 

The rest of the paper is structured as follows. First we specify our notation. Sec. [2] and [3] 
applies our expected-case analysis to estimating the real root separation bound, and to estimating 
the complexity of STURM solver. Sec. U] determines the expected number of real roots of random 
polynomial in the Bernstein basis and supports our bounds by experimental results. The paper 
concludes with a discussion of open questions. 

Notation. Og means bit complexity and the Og-notation means that we are ignoring 
logarithmic factors. For A = '^^^^ ai.X'' G 2[X], dg(A] denotes its degree. C{A) denotes an 
upper bound on the bitsize of the coefficients of A (including a bit for the sign). For a S Q, 
£ (a) > 1 is the maximum bitsize of the numerator and the denominator. A is the separation 
bound of A, that is the smallest distance between two (real or complex, depending on the 
context) roots of A. 

2 Subdivision-based solvers 

In order to make the presentation self-contained, we present in some detail the general scheme 
of the subdivision-based solvers. The pseudo-code of a such a solver is found in Alg. [TJ Our 
exposition follows closely [llj. 

The input is a square- free polynomial A G 2[x] and an interval Jq, that contains the real 
roots of A which we wish to isolate; usually it contains all the positive real roots of A. In what 
follows, except if explicitly stated otherwise, we consider only the roots (real and/or complex) 
of A with positive real part, since similar results could be obtained for roots with negative real 
part using the transformation x i— > —x. Our goal is to compute rational numbers between the 
real roots of A in Jq ■ 

The algorithm uses a stack Q that contains pairs of the form {f, J}. The semantics are that 
we want to isolate the real roots of f contained in interval J. PUSH(Q,{f, J}) inserts the pair {f, J} 
to the top of stack Q and POP(Q) returns the pair at the top of the stack and deletes it from 
Q. ADD(L, J) inserts J to the list L of the isolating intervals. 

There are 3 sub-algorithms with index SM, which have different specializations with respect 
to the subdivision method apphed, namely STURM, DESCARTES, or BERNSTEIN. Generally, 
INITIALIZATIONsM does the necessary pre-processing, COUNTsM(f)J) returns the number (or an 



Algorithm 1: SUBDIVISIONSolver(A, Jq) 


Input: Square- free A e 2M, 3o = [0, B] 

Output: A list of isolating intervals for the real roots of A in 


Jo 


1 1NITIALIZATIONsm(A, 3o) 




2 L^0, Q<-0, Q^ push(Q,{A, Jo}) 

3 while Q ^ do 

4 {f,a}<-pop(Q) 

5 V<— COUNTsM(f,J) 

6 switch V do 




7 




case V = continue 




8 




case V= 1 L<- add(L,J) 




9 




case V > 1 




10 
11 






{fL,JL},{fR,3R}^SPLIT™(f,3) 
Qf-PUSH(Q,{fL,aL}), Q<-PUSH{Q,{fR,jR}) 




12 RETURN L 





upper bound) of the real roots of f in 3, and SPLITsM(f, J) splits J to two equal subintervals and 
possibly modifies f . 

The complexity of the algorithm depends on the number of times the while-loop (Line [3] of 
Alg. [1]) is executed and on the cost of COUNTsM("f> J) and SPLITsM(f>3)- At every step, since 
we split the tested interval to two equal sub-intervals, we may assume that the bitsize of the 
endpoints is augmented by one bit. If we assume that the endpoints of Jq have bitsize T, then 
at step h, the bitsize of the endpoints of J C ^q is t + h. 

Let n be the number of roots with positive real part, and r the number of positive real roots, 
so r < n < d. Let the roots with positive real part, be aj = 9^(cXj) + iU(aj], where 1 < j < n 
and the index denotes an ordering on the real parts. Let A^ be the smallest distance between 
cci and another root of A, and Si = C (A^). Finally, let the separation bound, i.e. the smallest 
distance between two (possibly complex) roots of A be A and its bitsize be s = £ (A). 

2.1 Upper root bound 

Before applying a subdivision-based algorithm, we should compute a bound, B, on the (positive) 
roots. We will express this bound as a function of the bitsize of the separation bound and the 
degree of the polynomial. There are various bounds for the roots of a polynomial, e.g. [36^ 1161 ES] , 
and references therein. For our analysis we use the following bound [16] on the positive real 

„ V{k- 
2maXa^<omma^>o,k>i 



parts of the roots, B 



for which we have the estimation 



[miss] ccr < ^[ocn] <B< ^D\[(Xn). The bound can be computed in Oeld^T]. 

If we multiply the polynomial by x, then is a root. By definition of s, we have | log(|9^(ai) 
^(«j)l)l < s, for any i ^ j- Hence, we have the following inequalities 



IH(ai) - < 2 

m[oc2) - 9^(ai) < 2 



n0Cn-^) - noCn-2) < 2^ 

yi(an) - yt(cXn-l) < 2^ ( + ) 

IK(an) < n2^ 

Thus, we have B < j^in(an) < 5^rL2' < 16d^2' < d^2'*+\ Hence, we can deduce that 
/:(B) = 0(s + lgd). 



Lemma 2.1. Let A G Z[x], where dg(A) = d and C (A) = t. We can compute a bound, B, on 
the positive real parts of the roots of A, for which it holds B < d 2^^, and £ (B) = 0(5 + Ig d). 

Remark 2.2. In the worst case, the asymptotics of, more or less, all root hounds in the liter- 
ature, e.g. 136\ [T6l 126]/ . are same, since B < maxj, |ai,| < 2^ , and C (B) < T. However, it is very 
important in practice to have good initial bounds. Good initial estimations of the roots can 
speed up the implementation by 20% I22j. 

3 On expected complexity 

Expected complexity aims to capture and quantify the property for an algorithm to be fast for 
most inputs and slow for some rare instances of these inputs. Let E denote the set of inputs, and 
assume it is equipped with a probability measure (J.; then let c(I) denote the usual worst-case 
complexity of the considered algorithm for input I. By definition, the expected complexity is 
the integral J^c(I)|x(I). 

In our setting the set E depends on a parameter d (the degree of the input polynomial), 
and we are interested in the asymptotic expected complexity when d tends to infinity. Each 
E^i is equipped with a probability measure [iix (also called distribution) of the sequence of the 
(normalized) coefficients of the input polynomial and we consider the cases where there exists a 
limit distribution. 

3.1 Strategy and Independence 

A natural strategy is to decompose E^ into two subsets Ga and R^ (G stands for generic and R 
for rare), such that c(I) is small for I G G^ while 10.^(1) is very small for I G R^ and moreover 
the two partial integrals J^ c(I)|a.(i(I) and J^ c(I)|X(i(I) are balanced or at least both small. 

We face another difficulty. Classical properties and estimates in Probability theory are often 
expressed for a sequence of independent variables (i.i.d.) but most natural bijective transfor- 
mations performed in Computer Algebra do not respect independence. For instance, if X and 
Y are independent random variables, then U := X + Y and V := X — Y are not independent. In 
our setting, even if we consider a model of distribution of coefficients which assumes that they 
are i.i.d., then this does not imply that the roots are i.i.d. and we cannot apply usual tools or 
estimates. However, as we are interested in asymptotic behavior, for some models of distribu- 
tion of coefficients it happens that the limit distribution of the roots behave almost like a set of 
independent variables, i.e. they have very weak correlation. So we can invoke general classical 
estimates for our analysis. 

When this is not the case, a useful tool is the two-point, or multi-point, correlation function. 
They express the defect of independence between a set of random variables and classically serve, 
e.g., to compute standard deviations. 

Hereafter, we restrict ourselves to models of distribution of coefficients, hence induced dis- 
tribution of roots, for which the corresponding probability measures and correlation functions 
have already been studied. Hopefully these models will provide good approximations for the 
situations encountered in the many applications. 

3.2 S0(2) polynomials 

We consider the univariate polynomial A = ^^^q Qtx^, the coefficients of which are i.i.d. normals 
with mean zero and variances (^), where < 1 < d. Alternatively, we could consider A as 

A = ^1^0 V (i) '^i^^' where at are i.i.d. standard normals. These polynomials are considered 



by Edelman and Kostlan [S] to be "the more natural definition of a random polynomial". They 
are called S0(2) because the joint probability distribution of their zeros is S0(2) invariant, after 
homogenization. In |31) they are called binomial. Let p(t) = ,,VtJ] be the true density function, 
i.e. the expected number of real zeros per unit length at a point t G M. The expected number 
r of real roots of A is given by r = Jjp p(t)dt = yd [9j. Let aj be the real roots of A in their 
natural ordering, where 1 < j < r. 

We define the straightened zeros of A as 

Cj = Vidj) = vd arctan(aj)/7T, j = 1, . . . ,r, 

in bijective correspondence with the real roots ccj of the random polynomial, where V{t) = 
Jq p(u) du. Moreover, the ordering is preserved. The straightened zeros are uniformly dis- 
tributed on the circle of length 2vd |4j sec. 5]. This is a strong property and implies that 
the joint probability distribution density function of two, resp. ra, (distinct) straightened zeros 
coincides with their 2-point, resp. m-point, correlation function [1]. 

Proposition 3.1. [H Thm. 5.1] Following the previous notation, as d — ) oo the limit 2-point 
correlation of the straightened zeros is k(si , S2) — > 7t |si — S2I/4, when si — S2 — > 0. 

Let A(a) = mini<T^<T-{aT,+i — at} and A(C) = mini<i^<T-{Ci+i — Ci} be the separation bound 
of the real roots of A and the straightened zeros, respectively. We consider each straightened 
zero uniformly distributed on a straight-line interval of length 2vd- For two such zeros, we can 
consider one horizontal and one vertical such interval, defining a square, which represents their 
joint probability space. Since the real roots are naturally ordered, if two of them lie in a given 
infinitesimal interval, they must be consecutive. 

Let Z be a zone bounded above and below by a diagonal at vertical distance I from the main 
diagonal of the unit square. The probability Pr[A(C) < I] that there exist two zeros lying in 
a given interval of infinitesimal length I tends to the integral of k(si,S2) over the straightened 
zeros lying in Z, as d — > 00: 

Pr[A(0<l]^ J^k(si,S2)dsids2 

= 2j,'^dsi j:;+^k(si,S2)dS2 



Jo^dSlJ:;+^iSl-S2|dS2 = ^l^ 



where the first integral is over all straightened zeros, which lie in an interval of size 2vd- Notice 
that k(si, S2] is essentially the joint probability density function of two real roots. Using Markov's 
inequality, e.g. [28] we have Pr[A(C) > I] < E[A{Q]/\, so 

7T-2 ^fA 

E[A(C]] > I Pr[A(C) > I] = I - I Pr[A(C) < I] > I ^ ll 

This bounds the asymptotic expected separation conditioned on the hypothesis that it tends 
to zero, as d — > 00. If we choose I = l/(d'^T), where c > 1 is a (small) constant, which is in 
accordance with the assumption of I — ) 0, then E[A(C)] > jc:^ — t j3c"i/2 3 • 

E[A(0]=E[min{Ci+i-Ci}] = 

\/d l<i<T 

E[min{arctan(ai,) — arctan(ai+i)] = 

-- n l<i<r _ 

vd ,- . / at — ai+i Vi 1 7T 



E[min{arctan -— }] > — <=> 

7T i<i<T \1-|-aiai+iy d'^T 2d^'^"'/^T3 



b[min|arctan I )J > 



i<i<r \^ +0Ciai+iJ d'^+y^x 2d3cT-3- 

Function arctan is strongly monotone, and 1 +aiai.+i > 1 , for all i, except where oci is the largest 
negative root and at+i is the smallest positive root. But we can treat this case separately, since 
zero is an obvious separation point. 

E[ min {oCi - cxt+i }] > E[ min { f ""] ~ °''+' ) }] > 
1<i<r l<i<r \ I + a^ tti+i / 

7t 7T^ , 7t 7t^ 



> tan z-r^ -—3 — 7 > 



jJc+i/It- Zd^'^x^ d'^+y^T 2d3cT3' 
where the latter inequality follows from the series expansion tanx = x+x^/3+- • • for x G (0, 7r/2). 

Lemma 3.2. Let A € 2[x] of degree d, the coefficients of which are i.i.d. variables that follow 
a normal distribution with variances (.), then for the expected value of the separation bound 

of the real roots it holds E[A] > ,c+^/2 id^'^ ^ ' ^'^^ ^ constant c > 1, and E[s] = E[£ (A]] = 

0(lgd + lgT). 

3.3 Weyl polynomials 

We consider random polynomials, known as Weyl polynomials, which are of the form 



A = ^ aiX^/Vv., 



1=0 



where the coefficients Ui are independent standard normals. Alternatively, we could consider A 
as A = ^i^Q ciix\ where a^ are normals of mean zero and variance l/yl!. The density of the 
real roots of Weyl polynomials is 



^ ^ 1 / t2d(t2-d-l) t4d+2 

p(t) = -Wl + 



Tty et'r(Ti + i,t2) (et'r(Ti + i,t2))2' 

where V is the incomplete gamma function. The expected number of real roots is r = Jjp p(t)dt ~ 
- yd [31] , where the higher order terms of the number of real roots are not explicitly known up 
to now. 

The asymptotic density, for d — > oo, is 



pw= H ' r^x; (1) 




A useful observation is that the density of the real roots of the Weyl polynomials is similar 
to the density of the real eigenvalues of Ginibre random matrices, that is d x d matrices with 
elements Gaussian i.i.d. random variables [9l I31|. 

We consider only the real zeros of A that are inside the disc centered at the origin with radius 
yd since outside the disc there is only a constant number of them. In this case the density is 
represented by the first branch of ([1]). 

We work as in the case of the S0(2) polynomials. Now V{t) = Jq p(u)du = t/n. The 
straightened zeros, Ci, are given by 

Ci = V[(Xi) = (Xi/n, 

and they are uniformly distributed in [0, vd/7t] [31j . The joint probability distribution density 
function of two straightened zeros coincides with their 2-point correlation function. 



Proposition 3.3. |31| Under the previous notation, as d — ) oo the hmit 2-point correlation of 
the straightened zeros is w(si , S2) — > |si — S2l/(47r), when si — S2 — > 0. 

Working as in tlie case of tiie SO (2) polynomials, the probability Pr[A(C) < 1] that there 
exist two roots lying in a given interval of infinitesimal length I tends to the integral of w(si , S2) 
over the straightened zeros lying in Z, as d — > 00: 

Pr[A(0<l]= J2w(si,S2)dsids2 

= Jo^/"j:;!;w(s„S2)ds,ds2 = ^, 

and using Markov's inequality 

Pr[A(C) > I] < E[A(C)]/l ^^ E[A] > I- ^l^. 

If we choose I = 1/(d'^T), where c > 1 is a (small) constant, we get E[A(C)] > -^^ — lAic-yi 3 
and E[A(a]] > ^ - ^^^J-^/i^s - 

Lemma 3.4. Let A G 2[x] of degree d, the coefficients of which are i.i.d. variables that follow 
a normal distribution with variances 1/1!, then for the expected value of the separation bound 
of the real roots it holds E[A] > -^ - ^^^ic-yi^i and E[s] = E[£ (A)] = 0[lg d + Igx). 

3.4 The STURM solver 

Probably the first certified subdivision-based algorithm is the algorithm by Sturm, circa 1835, 
based on his theorem: In order to count the number of real roots of a polynomial in an interval, 
one evaluates a negative polynomial remainder sequence of the polynomial and its derivative 
over the left endpoint of the interval and counts the number of sign variations. We do the same 
for the right endpoint; the difference of sign variations is the number of real roots. 

We assume that the positive real roots are contained in [0, B] (Sec. 12. ip . If there are r of 
them, then we need to compute r — 1 separating points. The magnitude of the separation points 

subdivisions, performing 



is at most ^Aj, for 1 < j < r, and to compute each we need 



Ig^ 



binary search in the initial interval. Let T be the binary tree that corresponds to the execution 

of the algorithm and #(T) be the number of its nodes, or in other words the total number of 

subdivisions: 

T r 2B1 ^ 

#(T)=X Ig^ <2r + rlgB-^lgAj. (2) 

Using Lem. 12.1] we deduce that ^^(T) = 0{rs -|-rlg(d)]. 

The Sturm sequence should be evaluated over a rational number, the bitsize of which is at 
most the bitsize of the separation bound. Using fast algorithms |231I29| this cost is Oeld (t+s)); 
to derive the overall complexity we should multiply it by #(T]. Notice that for the evaluation 
we use the sequence of the quotients, which we computed in O^ld l] |231I29| . and not the whole 
Sturm sequence, which can be computed in O^ld^r], e.g. [7]. 

The previous discussion allows us to express the bit complexity of STURM not only as a func- 
tion of the degree and the bitsize, but also using the number of real roots and the (logarithm of) 
separation bound. This complexity is output sensitive, and is of independent interest, although 
it leads to a loose worst-case bound. 

Lemma 3.5. Let A G 2[x], dg(A) = d, C{A) = t and let s be the bitsize of its separation 
bound. Using STURM, we isolate the real roots of A with worst-case complexity Ob (rd (s +ts) ) , 
where r is the number of real roots. 



In the worst case s = 0{dr), and to derive the worst case complexity bound for STURM, 
Osld T ), we should also take into account that ds = 0[dT]. 

To derive the expected complexity we should consider two cases for the separation bound, 
that is, smaller or bigger than I = 1/(d'^T), where c > 1 is a small constant that shall be specified 
later. 

In the first case, that is A < I = l/(d'^T), the real roots are not well separated, so we rely 
on the worst case bound for isolating them, that is O^id'^T'^]. This occurs with probability 
Pr[A < I] = 0(vdl ) = 0( ,2c-i/2 2 ); by the computations of Sec. 13.21 and Sec J3.3l This 
probability is very small. 

For the second case, since A > l/(d'^T) we deduce s = 0(lgd + Igx). The complexity of 
isolating the real roots, following Lem. [33] is O^iTd r). The computations in Sec. l3.2l and Sec J3.3l 
suggest that this case occurs with probability Pr[A > I] = 1 —0[y/dl^] = 1 — 0( 2C-1/2 2 )1 which 
is close to one. 

The expected-case complexity bound of STURM is 

Ob (d - ^^H^) • rd^T+ ^J^ . dV) = Oelrd^T), 

for any c > 1, by using yd = 0{rr), which follows from the expected number of real roots. To 
avoid using this expected number, it suffices to set c > 2. 

Theorem 3.6. Let A e Z[x], where dg(A) = d, £ (A) = t. If A is either a S0(2) or a Weyl 
random polynomial, then the expected complexity of STURM solver is O-q [r d t) . 

In practice, the Sturm sequence is used and not the quotient sequence. The cost of the former 
is Oeld t) which dominates the bound of Th. 13.61 This explains the empirical observations that 
most of the execution time of STURM solver is spend on the construction of the Sturm sequence. 

4 Random Bernstein polynomials 

We compute the expected number of real roots of polynomials with random coefficients, repre- 
sented in the Bernstein basis. We start with some lemmata. 

Lemma 4.1. For k < n, non-negative integers, it holds 

j=0 ^ ' '' j=0 

Proof: We consider the RHS of the equality. For a specific j we expand the summand, and get 
terms of the form 



J^l^T^A^ kn-^gii^^^O<^<l^^_ 



There are kn + 1 such terms. Recall that e^^" = 1. Let \x. = Ak + v, where 1 < v < k — 1, 
< A < n, then 

"^^ \^kn-Ak-Vgi^(Ak+v) ^ 1^ TT-1^ \ kn-Ak-v ^i^^Ak gi^v ^ /^ "^^ \ kn-Ak-v gi^^v^ 



Ak + vy VAk + v/ VAk + v 

If we sum all these terms over j , we get 

^VAk + V VAk + V ^ 

9 





Figure 1. The transformation z : 



y+i 



since ^;=o e 



k-l „t^ 



0. 



Let \i = Ak. In this case, we have 



T^M kn-Ak gi^Ak 

Ak/ 



nk 
Ak 



^kn— Ak 



nk 

k(rL-A] 



Mn--\) 



Notice that < A < n. Summing up over ah A and all j, and multiplying by 1/k we get the 
LHS. D 



Lemma 4.2. For non-negative integers n, k, p, 



ill) 



n \p-i 



(^) 



1 



k(n-k) 



p-i 



Proof: The proof follows easily from Stirling's approximation n! ~ \/2nr\. [j)^- □ 

More accurate results could be obtained if the more precise expression \/2nr\. [j]^ e'2n+i < 
n! < \/2nn [j)^ efi^, is considered. 

4.1 The expected number of real roots 

We aim to count the real positive roots of a random polynomial in the Bernstein basis of degree 
d, i.e. 

p■.= ^bv(^^zH^-z]''-\ (3) 



where we assume that P(0)P(1 ] ^ 0, and {bi^} is an array of random real numbers, following the 
normal distribution, with "moderate" standard deviation, which shall be specified below. 

We introduce a suitable change of coordinates, z :=y/{y + ^), to transform a polynomial in 
the Bernstein basis into one in the monomial basis, by setting P = (1 +y)'^P(y). Now, P and P 
have the same number of real roots, and 

Even though the number of real roots does not change, their distribution over the real axis does, 
see Fig. [TJ In particular, we can now apply the techniques already used by Edelman, Kostlan, 
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and others for counting the number (and, eventually, the limit distribution) of real roots. Of 
course, by symmetry, the expected number of positive and negative real roots is equal. 
By Lem. 14.21 setting p = 2 and n = d we deduce: 



\ 



1 



7t V k(d-k)V Ul^ 




(4) 



It holds that -v/2/v7td < (i^)/'\/(2k) ~ V^ < 1 . To prove this, notice that Si^ is decreasing 
from 1 to d/2 and increasing from d/2 to d — 1 . Hence the lower bound is attained at k = d/2 
and the upper bound at k = 1 and k = d — 1 . 

Since S^ is small compared to (^) , it is reasonable to assume that omitting it will make only 
a negligible change in the asymptotic analysis. 

Let y = X , with x > 0. Now the problem at hand to count the positive real roots of 




k=o 



We need the following proposition 



Proposition 4.3. ^ Let v[t) = (fo(t), . . . , fn(t))''' be a vector of differentiahle functions and 
Co, . . . , Cn elements of a multivariate normal distribution with zero mean and covariance matrix 
C. The expected number of real zeros on an interval (or a measurable set) I of the equation 

Cofo(t) + ---+Cnfn(t)=0, is 



■ 1 

l7t 



|w'(t)||dt, w = w(t)/||w(t)| 



where w(t] = C^''^v(t). In logarithmic derivative notation, this is 



1 

7t 



92 



9x9ij 



log(v(x)TCv(x))k=,=tdt. 



For computing the integral in Prop. I4.3[ we shall use the logarithmic derivative notation. 
Following Prop. 1131 iidt) = y Cit)^'^^ and fii+i (t) = 0, Cii = ai and Cat+i = 0, where < i < d. 



and the variance is 1. Then, 



We consider function 



V X 



k=0 ^ ^ 



(xy) 



2k 



k=0 ^ ^ 



,2k 



,2k, 



By Lem. 14. 1[ for k = 2, we have f (z) = ^ ((1 + z 

d(z - 1 )2'i-i , f " (z) = d(2d - 1 ) (z + 1 )2'i-2 + d(2d - 1 ) (z - 1 )^'i. 

The following quantities are also relevant ff = ^d(z+ 1 )^'^^^ + dz(z 



z)^^) and so f'(z) = d(z+ Ij^'^-^ + 



m-'+\{z- 



l4d-l 



ld(2d - 1 ) (z + 1 )4d-2 + d(2d - 1 ) (z^ + 1 ) (z^ - 1 )2d-2 + 1 (2d - 1 ] (z - 1 )4d-2, and (f ]^ 



ff 

d^fz+"l)4d-2 + 2d2fz2 



ll^'i-i+d^fz-ll^i-^. 
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It holds that 



9x9y(logf(x,y)) 



f f + xyf f-xy(f 
f2 



'^2 



A 
f2' 



with 



A= d(z+l)4d-2(^(z+1) + ^(2d-1)z-zd) 

+d(z2 - 1 )2'i-2(z(z2 - 1 ) + (2d - 1 )z(z2 + 1 ) - 2d(z2 - 1 )z) 
-d(z - 1 f 'i-^l^lz - 1 ) + l(2d - 1 )z - zd) 



d ((Z + 1 )4d-2 + 4(2d - 1 )z(z2 - 1 )2d-2 - (z - 1 j^d-^) . 



If we let z = t , then 

A(t2) 



f(t 



TaT 



id((l+t2)td-2_|_4(2d-l)t2(t''-l]2d-2_(t2_i)4d-2) 
i((l+t2)2d+(l_t2)2d)2 



, l+(2d-l„-— 

TH 1 yi+t 



2t V(^-x^\^^-^ f^-t^V"-^ 



, l+t2 



. l+t2 



1 + 



We consider the substitutions t = tan j , tan 9 = t— p- , sin 6 = T—p- , cos 6 = t— ij , and 



^ — ^^,. Then 

A _ -,j 1 1+[2d-1)sin^9(cose)^^-^-(cos9) 

The expected number of positive real roots is given by 



I = irr^dt 



ttJo f(Fy 

1 pTT /J-T -^1 + (2d-l ) sin^ e (cos 9; 
^Jo V^O- l+(cos9^ 



d9 



IT" 



Performing the change 9 i— > tt— 6, we notice that I equals twice the integral between and 
7t/2. Hence, the expected number of positive real roots of P in (0,1) equals that in (l,oo). 
Hence, 

T — v^ r"/^ \/l+[2d-1]sin^9(cos9)^d-2_[cose)4d~ 
^~ 7t Jo l+(cos9)2d ^^■ 

Now we will bound the integral as d — > oo. Applying the triangular inequality and noticing 
that 1 + (cos9)4'i-2 < 1 and a + cos92'i > 1, we get 



I < 



jo^W^de + J7^ V2d^^sin9(cos9; 



d-l 



d9 



7T d 



y/Id 
2 



- /- 

TtV d 



2 "•" 7t- 



For a lower bound, we neglect the positive term (2d — 1) sin 9 (cos 9)'^'^ ^, and notice that 
^1 + (cos9)4d-2 > 1 + (cos9)2'l-i > 1 + (cos9]2'i-^ = (1 + (cos9)'l-l)(1 - (cos9)'i-M, and 
> 1. 



1+(cos9)''~' 
l+(cos9)2d 

Lemma 4.4. 



Win] 



'7t/2 9 

(cos9]"-d9 < 





1 



Proof: We need the following inequality [6J on Wallis' cosine formula: 

1 1 • 3 • 5 • • • (2k - 1 1 1 



1 • 3 • 5 • • • (2k - 1 ; 
< !^ < 



V^7r(k + 47T-1 - 1 ) " 2 • 4 • 6 • • • (2k) " ^7t(k+l/4) 
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If n :s even then Win) = f M!l = fS^ < 7 =;T) ^ ^7^ " 
If n :s odd then W[n) = ^^^ = '■','^tiSt^' < V^^^^ + ^^^ " ^ ) " ^ < ^tJtT' ° 
Using the lemma, Jq (cos9)'^^^ d9 < ^-7^, so: 



2 ^ '-4^;^ 2^ 



Hence I = -^^ ±0(1) and we can state the fohowing: 



Theorem 4.5. The expected number of real roots of a random polynomial P = ^j^^q "^l^V (ik) 
where a^ are standard normals i.i.d. random variables, is v2dib 0(1). 



2d\ 7k 



By employing (|4]) and considering \/Sk as part of the deviation, we have the fohowing: 

Corollary 4.6. The expected number of real roots of a random polynomial in the Bernstein 
basis, Eq. (0), the coefficients of which are normal i.i.d. random variables with mean and 

variance 1/Sk = Vxhn^^, is ^/2d±0[^). 



7rk(d-k) ■ 

In Table [1] we present the results of experiments with polynomials in the Bernstein basis 
(see Eq. ([3])), of degree < 1000, the coefficients of which are i.i.d. random variables following 
the standard normal distribution, that is mean zero and variance 1. For each degree we tested 
100 polynomials. The first column is the degree, while the second is the expected number of 
real roots predicted by Cor. 14.61 which assumes variance 1/Si^. The third column is the average 
number of real roots computed. Our experiments support the following conjecture: 

Conjecture 4.7. The expected number of real roots of a random polynomial in the Bernstein 
basis, Eq. 0, the coefRcients of which are standard normal i.i.d. random variables, that is with 
mean and variance 1 , is \J2A it O f 1 1 . 



Columns 4-7 of Tab. [T] corresponds to the average number of real roots in the intervals 
(— 00,— 1), (—1,0), (0,1) and (l,oo), respectively. For these experiments we took random poly- 
nomials in the monomial basis and converted them to the Bernstein basis. The roots of a 
random polynomial in the monomial basis, under the assumptions of [T7], concentrate around 
the unit circle. The symmetry of the density suggests that each of the intervals (—1/(1 — p), — 1], 
(— 1,— l-|-p), (1— p, 1), and (1,1/(1 — p)), contains on the average 1/4 of the real roots (Fig. [H 
left). If we apply the transformation x 1— > x/(x -|- 1 ) (Fig. [H right) to transform the polynomial 
to the Bernstein basis, then 3/4 of the real roots are positive, 1 /2 of them are in (0, 1 ) and 1 /4 
in (1, 00). We refer to the last columns of Tab. [T]for experimental evidences of this. 

As far as the distribution of the real roots in (0, 1 ) is concerned, if we denote them by tt, 
then arccos(2ti — 1], behaves as the uniform distribution in (0, 7t]. In Fig. [21 we present the 
probability-probability plot, (using the ProbabilityPlot command of MAPLe) of this function 
of real roots of random polynomials in Bernstein basis, of degree 1000 (light grey line), against 
the theoretical uniform distribution (black line) in (0, 7t). We observe that the lines almost 
match. For reasons of space, we postpone the discussion about the distribution of the roots. 
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Figure 2. Left: Function arccos(2t — 1) of real roots in (0, 1) 
polynomials in the Bernstein basis for d £ {5, 10, 15}. 



against uniform distribution in (0, 7t). Right: Density of 



d 


yjd 


( — oo, oo) 


(-00,-1) 


(-1.0) 


(0,1) 


(l,oo) 


100 


14.142 


13.640 


0.760 


2.740 


6.530 


3.610 


150 


17.321 


16.540 


0.890 


3.260 


8.090 


4.300 


200 


20.000 


19.740 


1.100 


3.780 


9.740 


5.120 


250 


22.361 


21.400 


1.350 


3.970 


10.610 


5.470 


300 


24.495 


24.320 


1.270 


4.760 


12.300 


5.990 


350 


26.458 


26.540 


1.620 


5.100 


13.400 


6.420 


400 


28.284 


27.980 


1.490 


5.430 


14.080 


6.980 


450 


30.000 


29.460 


1.620 


5.890 


14.970 


6.980 


500 


31.623 


31.200 


1.830 


5.960 


15.620 


7.790 


550 


33.166 


32.740 


1.770 


6.360 


16.290 


8.320 


600 


34.641 


34.300 


1.850 


6.570 


17.270 


8.610 


650 


36.056 


35.480 


2.050 


6.840 


17.240 


9.350 


700 


37.417 


37.200 


2.160 


7.510 


18.650 


8.880 


750 


38.730 


38.180 


2.190 


7.300 


19.360 


9.330 


800 


40.000 


39.160 


2.220 


7.830 


19.490 


9.620 


850 


41.231 


40.420 


2.130 


8.010 


20.320 


9.960 


900 


42.426 


41.780 


2.390 


8.070 


20.530 


10.790 


950 


43.589 


42.680 


2.200 


8.330 


21.570 


10.580 


1000 


44.721 


43.540 


2.400 


8.610 


21.770 


10.760 



Table 1. Experiments with random polynomial in the Bernstein basis. 

5 Conclusions and future work 

Our results explain why the solvers are fast in general, since typically there are few real roots 
and in general the separation bound is good enough. This agrees with the fact that in most 
cases the practical complexity of the STURM solver is dominated by the computation of the 
sequence and not by the evaluation. Our current work extends the first part of this paper to 
Kac polynomials, and to solvers DESCARTES and BERNSTEIN. 

The main issue with the Kac polynomials is that there is a discontinuity at ±1 when d — > oo. 
To be more precise, the fact that there are few roots even near ±1 , where they are concentrated 
asymptotically, is balanced by the fact the 2-point correlation, k(si , si), between two consecutive 
roots is a complicated function of |si — S2I, Si and d and (in opposition with the two other 
distributions we studied) its limit when d tends to infinity is not equivalent to a simple function 
of |si — Sil- This is an interesting problem which deserves to be studied and investigate further. 

An interesting question is whether we can design a randomized exact algorithm based on 
the properties of random polynomials. Lastly, we wish to extend our study to polynomials with 
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inexact coefficients. 
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